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ABSTRACT 



Aims. We present an unbiased orbit solution and mass determination of the components of the eclipsing binary PG 1336-018 as a 
critical test for the formation scenarios of subdwarf B stars. 

Methods. We obtained high-resolution time series VLT/UVES spectra and high-speed multicolour VLT/ULTRACAM photometric 
observations of PG 1 336-0 1 8, a rapidly pulsating subdwarf B star in a short period eclipsing binary. 

Results. Combining the radial velocity curve obtained from the VLT/UVES spectra with the VLT/ULTRACAM multicolour 
lightcurves, we determined numerical orbital solutions for this eclipsing binary. Due to the large number of free parameters and 
their strong correlations, no unique solution could be found, only families of solutions. We present three solutions of equal statistical 
significance, two of which are compatible with the primary having gone through a core He-flash and a common-envelope phase de- 
scribed by the a-formalism. These two models have an sdB primary of 0.466 M and 0.389 M Q , respectively. Finally, we report the 
detection of the Rossiter-McLaughlin effect for PG 1336-018. 
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1. Introduction 

The subdwarf B (sdB) stars are generally acknowledged to 
be core helium burning stars with a canonical mass of 
approximately 0.5 M G . Their thin, inert hydrogen envelope 
(M em < 0.02 M ) places them on the hot extension of the 
Horizontal Branch (HB), the so-called Extreme Horizontal 
Branch (EHB). Since the hydrogen envelope is too thin to sustain 
nuclear burning, these stars will not go through the Asymptotic 
Giant Branch and Planetary Nebula phases. Instead, when their 
core helium has run out, they will enter a He-shell burning phase, 
where they expand and heat up, making them appear as sdO 
stars before they evolve directly onto the white dwarf cooling 
sequence. Even though the models describing the future evo- 
lution of the sdB s tars are generally accepted (e.g. those of 
iDorman et alJ Tl993). the current evolutionary state of the sdB 
stars is still poorly understood. The fact that sdB stars must have 
lost almost all of their hydrogen layer at exactly the same time 
when the helium core has attained the minimum mass required 
for the helium flash to occur, makes them enigmatic from an evo- 
lutionary point of view. To loose such an amount of mass, they 
must suffer considerable mass loss during the red giant branch 
(RGB) phase, and most probably also during the helium core 
flash. The most fundamental missing piece to our understanding 
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of the evolution of the sdB stars, apart from the physics during 
the h elium core flash, is the natur e and physics behind this mass 
loss dFusi-Pecci & Renzinilll976J) . 

In recent years it has been discovered tha t a sign ificant frac- 
tion of sdBs are in binaries. iMaxted et al.l d200ll) found that 
about t wo-thirds of the sdB star s in the field are members of bi- 
naries. Napi wotzki et al.l (|2004|) found a binary fraction of 40% 
among stars i n the SPY (Supernova type la Progenitor) survey 
sample, while Mor ales-Rueda et all (|2006) found 48% in a sam- 
ple from the Edinburgh-Cape (EC) survey. Many of the binary 
sdBs are found to be in short period systems with periods from 
a few hours to several days, with companions bei ng either white 
dwarfs or M-dwarfs dMorales-Rueda et alj |2003). The peculiar 
frequency of binarity has been an important constraint on evolu- 
tionary population synthesis theory, and has led to the acknowl- 
edgment that the binarity has to play a key role in the formation 
channels f or sdB stars. There are several binary mechanisms pro- 
posed by (lHan et ai1l2002l 120031, and references therein) as for- 
mation channels for sdB stars : 



1. common envelope ejection, leading to short-period binaries 
with periods between 0.1 and 10 days and an sdB star with 
a very thin hydrogen envelope, and with a mass distribution 
that peaks sharply at 0.46 M Q . Depending on the secondary, 
a main-sequence star or a white dwarf, the subchannels are 
called the first CE ejection channel and the second CE ejec- 
tion channel, respectively, 
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2. stable Roche lobe overflow, resulting in similar masses as in 
1 . but with a rather thick hydrogen-rich envelope and longer 
orbital periods between 10 and 100 days, 

3. double helium white dwarf mergers giving rise to single sdB 
stars with a wider distribution of masses. 

Detailed investigation of sdB binaries is crucial in order to 
determine their masses for comparison with the theoretically 
proposed evolutionary channels. New momentum in the efforts 
to resolve the evolutionary paths of sdB stars ca me a decade ago, 
after the discovery that some of them pulsate ( Kil kenny et al.l 
1 19971) . This has opened up a new window into their interiors via 
the techniques of asteroseismology and stimulated a burst of re- 
search. Extensive search campaigns have revealed two classes of 
pulsating sdB stars known as short period sdB variables (sdBV 
or V361 Hya stars, formerly EC 14026 stars, after the prototype) 
and long p eriod sdB variable s known as PG 1716 stars (or lps- 
dB V stars, iGreen et al.ll2003l) . 

The sdBV stars, discovered bv lKilkennvetalld 19971) and in- 
dependently theoretically predicted by ICharpinet et aL (Q996), 
are low amplitude multimode pulsators with typical periods 
ranging between 100-250 s. Their pulsation amplitudes are gen- 
erally of the order of a few hundredths of a magnitude. The short 
periods, being of the order of and shorter than the radial funda- 
mental mode for these stars, sugge st that the observed m odes 
are low-order, low-degree /j-modes (ICharpinet et al.l l2000). The 
39 known sdBV stars occupy a region in the r e ff-logg plane 
with effective temperatures between 28 000 K and 36 000 K and 
surface gravities (logg) between 5.2 and 6.2. 

The detailed asteroseismological modelling of sdBV stars 
is hampered by the fact that there are too few pulsational fre- 
quencies to fit t hose predicted from non-rotating or rigidly 
rotating models dBrassard et all 120011: ICharpinet etaTI 120051: 
iRandall et al.l 12005 ). However, the observed frequency spectra 
are too dense to be accounted for by only low-degree (£ < 2) 
modes. In order to have a unique asteroseismological model we 
need to have accurate pulsation frequencies and an unambiguous 
identification of the modes of oscillation (spherical wavenum- 
bers £ and m). Thanks to multisite campaigns by the WElQ de- 
voted to resolving the frequency spectrum of sdBV stars in the 
last decade, we do have extensive and reliable frequency lists 
for several sdBVs. The problem lies in the second requirement 
mentioned above, the unambiguous mode identification. There 
are only t wo ways this can be achieved: through line profile 
variations (jAe rts & Ever 2000) or the amplitude ratio method 
dDupret et alJl2003URandall et al.ll2005l) . 

As sdBV stars are quite faint (the brightest one is m% = 1 1.8) 
and their periods are very short, the line profile variation method 
poses a real challenge considering the low S/N that accompa- 
nies any high-resolution time-resolved spectroscopy, even with 
the biggest telescopes available. Hence, the line profile varia- 
tion method has not yet been reliably applied to any sdBV star. 
The amplitude ratio method is not problem free either. Due to 
the very low pulsational amplitudes, the photometric errors are 
usually too large for unambiguous identification of the spherical 
degree i of th e modes, especially to distinguish between t- 0, 1 
and 2 modes dJefferv et al.ll2005l) . 

Among the binary sdB stars, four eclipsing sdB sys- 
tems have been discovered that all show a deep and 
strong reflection effect, with very short orbital periods in 
the rather narrow range of 130-170minutes. Such short or- 
bital periods imply that they must have evolved through 
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binary mass transfer and common envelo pe evolution. Out 
of these four systems, namely HW Vir dWood et al.l [1993; 
iMenzies & Marangll 19861) . NY Vir (iKilkennv et al.ll 19981) (h ere- 
after PG 1336- 018), HS 0705+6700 dDrechsel et alJ l2001h and 
HS 2231+2441 d0stensen et al.ll2007l) . only one system contains 
a rapidly pulsating sdB star as a primary: PG 1336-018. As 
such, this system provides a natural laboratory for detailed evo- 
lutionary and asteroseismic analyses, which is the purpose of our 
project. 

PG 1336-0 18 was classified as an sdB star in the Palomar- 
Green survey dGreen et al ] 119861) and shown to be a close 
ecl ipsing binary with sh ort-period multimode light variations 
by IKilkennv et al.l dl998l) . Assu ming the primary mass to be the 
canonical sdB mass of 0.5 M , IKilkennv et al.1 (1 19981) find that 
the secondary must be a mid-M dwarf with a mass of about 
0.15 M Q . Soon after its discovery, PG 1336-018 was a target of 
two Whole Earth Telescope (WET) campaigns, Xcov 17 in April 
1999 dKilkennv et al J 12003b and Xcov 21 in April 2001. These 
white light data resolved more than 20 frequencies in the tem- 
poral spectrum dKilkennv et al.ll2003b in the range from 5000 to 
8000 //Hz. Even though the frequency content of the star is thus 
known very precisely, an adequate asteroseismic model is still 
lacking mainly due to the lack of an unambiguous mode identi- 
fication. The colour behaviour is needed for photometric mode 
identification to identify the spherical degree I of the modes and 
to discriminate between the numerous possible seismic models. 
To further reduce the allowable seismic model space we need 
to examine line profile variations due to the pulsations in order 
to disentangle the azimuthal wavenumber m. Only with the ac- 
curate pulsation frequencies and an unambiguous mode identi- 
fication can the asteroseismology provide the accurate mass es- 
timate needed for confrontation with those predicted from the 
formation scenarios for sdB stars. 

PG 1336-018, being the only rapidly pulsating sdB star in 
an eclipsing binary, is the only star with enough potential to 
confront the proposed evolutionary scenarios, as the eclipses 
help constrain the inclination and radii. Therefore we study 
PG 1336-018, this time armed with new multicolour photomet- 
ric and spectroscopic VLT data. In this paper we present the new 
data and the orbital solution. This is the first step toward our ulti- 
mate goal, an accurate mass determination of PG 1336-018 and 
a critical test of current stellar evolution theory. 



2. Observations and data reduction 

2.1. Photometry 

PG 1336-018 (a 2 ooo = 13:38:48.2, foooo =-02:01:49.0, m v = 
13.4) was observed on the night of May 18/19 2005 using the 
ULTRACAM camera attached to the ESO VLT UT3 (Melipal) at 
Paranal Observatory, Chile. ULTRACAM is a high-speed three- 
channel CCD camera specifically de signed for fast photometry 
programmes dDhillon & Marshll200lb . We gathered two full or- 
bital cycles, about 5 h, of PG 1336-018 si multaneously in thre e 
filters f, g' and u' of the SDSS system dFukugita et al.l lT996). 
The seeing (around 0.9 arcsec) was variable during the night 
and getting worse toward the end of the run. The exposure time 
was 0.5 s in the beginning of the run, but due to poorer seeing 
was increased to 1 s to improve the S/N. This did not deteriorate 
our temporal resolution signi ficantly, since the sho rtest period 
found in PG 1336-018 is 97 s dKilkennv et alJl2003l) . To achieve 
1 second time resolution, it was necessary to define 2 windows 
on each of the 3 ULTRACAM chips. One window was placed 
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Fig. 1. ULTRACAM/VLT r (upper), g' (middle) and m' (bottom) lightcurves of the eclipsing sdBV star PG 1336-018 from 2005 
May 18/19. The insets show enlarged sections of the two primary eclipses, where pulsations are clearly visible. The differences 
between the two consecutive primary eclipses, apart from the noise, are due to the beating of the modes and different phases covered 
during the eclipse. The shape of the li lightcurve is discussed in the text. The ordinate is the differential magnitude, and the abscissa 
is Fractional Julian Date. 



around PG1338-018, and another on a nearby comparison star. 
The dead-time of the observation was 24 milliseconds. 

All data frames were reduced using the UL TRACAM 
pipeline reduction software (Dhillon & Marsh 2001). Care was 
taken to select the most optimal choices offered in the reduc- 
tion software. The 'normal' extraction method with the 'vari- 
able' aperture sizes, as they track local changes in the seeing 
disk, gave the best results. Several apertures were tried out and 
an aperture of 1.7 times the FWHM gave the highest S/N for 
r' and g' band. The star counts were divided by the compar- 
ison star counts and converted to obtain a differential magni- 
tude (V-C) in each filter. As both the target and the comparison 
star were in the same field, differential photometry accounted 
well for the variations in the sky transparency and extinction in 
r' and g' band. Unfortunately, the only comparison star within 
ULTRACAM's 2.6 arcminute field of view on the VLT is very 
faint in the blue, resulting in poorer differential photometry in 
the m' compared to the f and g' band. Therefore, a wider aper- 
ture had to be used for the li band. Due to the faintness of the 
comparison star in w', its g' band lightcurve was used to make 
the differential u' lightcurve. This gave a satisfactory result in 
the sense that both the pulsations and the eclipses were recov- 
ered, but it introduced an unreliable slope in the first part of the 
m' lightcurve (see Fig. [T). Therefore, we did not rely on the m' 
lightcurve for the orbital analysis. However, we did use the sec- 
ond part of the m' lightcurve to cross-check our results, as well 
as for the frequency analysis (see Sect. I4.21 i. 

The times in the data frames were converted to JD and 
barycentrically corrected. Differential (V-C) lightcurves for r', 
g' and m' were constructed from a set of more than 80000 sci- 
ence frames. The f , g' and m' lightcurves are plotted in Fig.Q] 
where we can see a clear sign of the pulsations of the primary 
component in all the phases of the binary orbit, even during the 
primary eclipse. A strong reflection-like effect (0.2 magnitudes 
in g' and 0.25 magnitudes in r') is evident. This effect, charac- 
teristic of all binary systems containing an sdB star and a cool 
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Fig. 2. A typical single UVES/VLT spectrum of PG 1336-018 
from our VLT run on 2005 April 28 (top) and the coadded 
spectrum (bottom), produced by combining all the 399 avail- 
able spectra after shifting according to the orbital radial veloc- 
ity solution. The Balmer lines are indicated together with the 
helium lines used for the determination of physical parameters. 
Discontinuities due to imperfect merging of spectral orders only 
become evident in the high-S/N combined spectrum. 



M-dwarf companion in rotationally locked orbit, is due to the 
high contrast in the temperatures between the heated and un- 
heated hemispheres of the M-dwarf. 



2.2. Spectroscopy 

Even though PG 1336-018 was a target of several photomet- 
ric campaigns, its faintness relative to the rapid oscillations 
has prevented any reasonably good time-resolved spectroscopy. 
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Fig. 3. A sample fitting of two Gaussians to the observed H y line 
(the same spectrum as the one shown in Fig. [2j using molly. 
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Fig. 4. The radial velocity measurements (average of the H E , H,5, 
H y and lines) of all the individual UVES/VLT spectra. The 
best fit orbit solution from PHOEBE is also shown. 



The short pulsation periods require very short integration times. 
There were two attempt s o far with the aim o f detecting the pul- 
sational radial velocities dWoolf etalJ 120031) and identification 
of the pulsa tion modes from the wavelength dependency of the 
amplitudes dDreizler et alJl2000h . both with a null result. 

A time-series of 399 high resolution spectra were taken over 
a period of ~ 9 h, covering about 3.7 full orbits, on the night of 
April 28, 2005 using the Ultraviolet Visual Echelle Spectrograph 
(UVES) on the VLT UT2 (Kueyen) at the Paranal Observatory, 
Chile. Only the blue arm was used, with wavelength coverage 
from 3900 to 5000 A, and the slit width of 1 arcsec at a resolution 
of 46 890. Each spectrum was integrated for 45 s which, with the 
ultra fast read-out of about 23 s we used, gave a time resolution 
of 68 s. Dome flat-fields and bias calibration frames were taken 
at the beginning and at the end of the night, and ThAr exposures 
were taken before and after the run. 

Due to the very low signal we got for such a short expo- 
sure and the ultra fast read-out mode used, the UVES reduction 
pipeline did not give satisfactory results. Therefore, we devel- 
oped a non-standard reduction method, using the ESO-MIDAS 
package. This provided a factor of ~ 2 increase in the S/N ratio of 
the reduced spectra, compared to those produced by the pipeline. 
The bias calibration frames had an offset between the upper and 
the lower part, due to the ultra fast read-out mode used. After 
careful examination of each bias frame, we proceeded as fol- 
lows. First we examined the interorder space of each science 
frame (by taking the median of the box) to determine these off- 
sets which were then subtracted from the science frames. Then 
the science frame was corrected for cosmic rays, extracted and 
background corrected (which was smoothed to reduce the noise). 
Since, in our case, the sky background contributes most to the 
noise, we used optimal extraction which gave better S/N, as sug- 
gested bv lMukaa dl990h . Then the science frames were fiat-field 
corrected, wavelength calibrated and, finally, the orders were 
merged. Since the spectra were oversampled we have rebinned 
them in an optimal way such that the S/N increased without com- 
promising the resolution. Finally, the science frames were nor- 
malized. 

A typical individual spectrum of PG 1336-018 is shown in 
the top panel of Fig. [2] The bottom panel of Fig. [2] shows the 
coadded orbit-corrected spectrum (see Sect. 14.11 ). Despite our 
extensive effort to achieve the optimal reduction scheme, the ex- 
traction and merging of the orders is not perfect. This is due 



to the fact that the Echelle order discontinuities do not behave 
'consistently' under a low signal. This leads to some jumps and 
wiggles seen in the continuum of the coadded spectrum and par- 
ticularly in the red wing of H y . For this reason we did not make 
use of this line in the merged spectrum for the spectroscopic pa- 
rameter determination discussed below. 

In the blue wavelength range covered by our data no sign 
of any spectral feature fr om the cool comp anion can be seen, 
confirming the results of IWoolf et alj d2003l) . Due to the large 
difference in effective temperatures (about a factor of 10, see 
Sect. H]i the hot sdBV star dominates the spectrum even in the 
primary eclipse. 

3. RV determination 

Our spectra allow us to produce a radial velocity (RV) curve, 
with an excellent phase coverage, from which we can indepen- 
dently determine the orbital period (P) and semi-amplitude (K\) 
of this eclipsing binary. As we are dealing with a low S/N, we de- 
termined RVs from the spectra trying out several different meth- 
ods. The best results were obtained by using molly- a software 
package, which fits two Gaussian profiles to the Balmer line pro- 
files [jf This allows good treatment of both the broad wings and 
the sharper core at the same time. This gave better results than 
any of the other methods we have tried. 

We have measured the RVs of the highest S/N lines in the 
spectrum, namely H e , H,5, H r and Hp, using this package. A sam- 
ple fit is shown in Fig.[3]for an individual spectrum. The FWHM 
of the two Gaussian fits, as well as their heights, were treated 
as a free parameter at first, but were kept fixed once the best fit 
values were found. We checked carefully if the RV from the H y 
line deviated from the one of the other Balmer lines, due to the 
discontinuity in its red wing. This turned out not to be the case 
(see also Fig. [3) so we kept the H r RV values in our analysis. 

Finally, the average of each RV measurement, using H e , H^, 
H y and Fig lines, was determined. These radial velocity values 
for each of the 399 individual spectra (with the errors), are shown 
in Fig. |4] together with the best fit orbital solution (see Sect.|4|. 

To perform an independent determination of the orbit from 
our spectroscopic data, and to verify the photometric ephemeris, 

2 http://deneb.astro.warwick.ac.uk/phsaap/software/molly/html/IN- 
DEX.html 
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the measured RVs (after barycentric correction of the velocities 
and the mid-exposure times) were subje cted to a periodogram 
analysis. A sinusoidal fit using Period04 dLenz & Breger 2004) 
gives the frequency 1 14.25 + 0.1 fiHz and the semi-amplitude 
78.6 ±0.6 km/s which is, considering our poor frequency reso- 
lution of about 30 /vHz, in a good agre ement with the orbital pe- 
riod P= 0.101015999 d calculated by IKilkennv et al.l d2000t) as 
well as with the values derived from PHOEBE in Section [4] The 
semi-amplitude of the velocity var iation is in goo d agree ment 
with the 78 ± 3 km/s estimated by IKilkennv et all ([1998) (see 
their Table 4) even though they reported the semi-amplitudes of 
all of their observations (see their Table 3) to range from 47 + 4 
to 79 + 4 km/s. The semi -amplitude we obta ined is somewhat 
larger than estimated by IWoolf et af] d2003lh 64 ± 1 km/s, but 
their data cover only 1 .4 orbits and contain a gap which prob- 
ably resulted in an underestimated value. 

As our data set suffers from a baseline too short for reliable 
ephemeris determin ation, we adopted the ephemeris obtained by 
IKilkennv et all d2000t) (see TablefTT). 

Since the system is single-lined and the orbit is assumed to 
be circular, the analysis of the RV curve is straightforward. The 
mass function calculated from the semi-amplitude and the period 
gives: 

/(M) = 0.0051 +0.0001 M . 

4. Orbital parameters 

In order to investigate the pulsational properties of 
PG 1336-018, the subject of a follow-up paper, the orbital 
variations due to the binarity must be removed from the ob- 
served lightcurve. However, in order to find the best orbital 
solution for this eclipsing binary system, the pulsations of the 
sdB primary must be removed as well. This is a non-trivial 
coupled problem. The determination of the orbital parame- 
ters of this system required to understand and evaluate the 
temporal spectrum of the primary sdB pulsator. In order to 
achieve this, we followed an iterative procedure, using all the 
information about the target we have. Once we find a reliable 
orbital solution, we subtract it from the lightcurves. Then we 
use the orbit subtracted lightcurves to extract the pulsation 
frequencies present in our data. We prewhiten the original 
observed lightcurves with these frequencies. The prewhitened 
lightcurves are then used as input to find the second iteration 
orbital solution. 

4.1. Fundamental parameters 

Our high resolution VLT/UVES spectra allow us to improve th e 
spectroscopic parameters determined by IKilkennv et alj (Q998). 
Using our RV solution (see Fig. 0J, we shifted the spectra and 
added them together to improve the S/N. The coadded orbit- 
subtracted spectrum is shown in the bottom panel of Fig. [2] 

Fo r the model fittin g procedure, we used the LTE mod- 
els of iHeber et alj d2000l) . The model spectra were convolved 
with a Gaussian instrumental profile of 0.25 A and rotationally 
broadened (assuming tidally locked rotation) with a vsinj of 
74.2 km/s. This produces a model spectrum with line cores that 
reproduce the observed spectrum excellently for all lines that 
are unaffected by Echelle order discontinuities. Unfortunately, 
while the fit to the cores is good, the wings are not well fitted. 
Our best simultaneous fit for effective temperature, gravity and 
helium abundance yields: 

r eff = 3 1 300 + 250 K 




log He/H = -2.929±0.009 



-50 50 

AX (A) 

Fig. 5. Our spectroscopic model fit to the mean spectrum in 
Fig. |2 The best fit model spectrum has been plotted on top of 
the observed spectrum as a smooth curve. Note that the Hy line 
was kept out of the fit due to its proximity to an echelle order 
discontinuity. 



logg = 5.60 ± 0.05 dex 
logy = -2.93 ± 0.05 dex 

The quoted errors are about five times larger than the formal fit- 
ting errors reported in Fig. [5] Although such 5<x errors would 
normally be quite conservative considering the resolution and 
signal of the combined spectrum, there are obvious problems. 
The effects of errors due to the Echelle extraction problems de- 
scribed earlier are hard to quantify. The effective temperature is 
well constrained by the depth of the high order Balmer lines, and 
the helium abundance is determined by the depth of the narrow 
He i lines (marked in Figs.|2]and|5]i, which are not much affected 
by the Echelle extraction problems. However, since the Echelle 
order discontinuities strongly affect the wings of the lines, which 
are essential for the gravity determination, we cannot exclude a 
large error on log g. For this reason, we will only use the effective 
temperature determination as a constraint for our orbital fitting 
procedure, and not log g. Indeed, as we will see later, such a low 
log g is inconsistent with any realistic mass-radius relationship 
that can be derived from the orbit by at least 0.15 dex. In order 
to rule out other causes for the inconsistent log ^determination 
from the average spectrum, we tried to fit it using NLTE at- 
mosphere models, enhanced metallicity models, or changing the 
assumed rotational velocity broadening. All these attempts pro- 
duced negligible changes to the derived parameters listed above. 
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Fig. 6. The ULTRACAM/VLT g' lightcurve together with the synthetic orbit solution. The middle panel shows the residuals of the 
orbit subtraction. Pulsations during the eclipses are now clearly visible, and we can see that the amplitude is smaller during the 
primary eclipse than during the secondary as only the part of the surface is visible. The bottom panel shows the residuals after 
prewhitening with the four strongest oscillation modes. 



4.2. Binarity and pulsation 

Numerical orb it solutions were in vestigated using the PHOEBE 
package tool (IPrsa & Zwitterl 120051) which incorporates the as- 
pects of the Wilson-Devinney (WD) code dWilson & DevinnevI 
1197 lb . The WD approach uses differential correction (DC) as 
the minimization method, which is in essence a linearised least 
squares method. The code was used in the mode for detached bi- 
naries with no constraints on the stellar potentials. No third light 
or spots were included. 

The ULTRACAM/VLT g and r' lightcurves and the RV 
measurements obtained from the UVES/VLT spectra were 
solved simultaneously to yield a consistent model fit. As PHOEBE 
is limited by the number of points (currently the limit is 9000 
points) we had to phase bin our ULTRACAM/VLT lightcurves 
into 4000 data points per lightcurve. 

The major problem in finding the orbital solution of any bi- 
nary system is not only the fact that there are many free param- 
eters (12 + 5n, where n is the number of lightcurves in different 
filters), but also that the parameters are correlated. Some of these 
correlations are severe, especially between the mass ratio q and 
the potential of the secondary star Q2 (see the discussion be- 
low in Sect. 14.3b . Hence, one is left with several formal families 
of solutions within the parameter space. We must then confine 
the range of possible solutions by reducing the number of free 
parameters. The only safe way to do this is by considering the 
boundary conditions set by the data themselves and by sound 
theoretical considerations. 

The parameters that were assumed and kept fixed in our anal- 
ysis were ?o, P, T e s of the primary, gravity darkening coefficients 



Table 1. Fixed parameters in the search for the orbital solution 
of PG 1336-018. 



Parameter 


Value 


to 


2450223.36134 d a 


P 


0. 101015999 d a 


T e m 


3 1300 K 


T e m 


3000 K b 


gi 


1.0 


g2 


0.32 


Ai 


1.0 


M (gl 


0.217 


x-i (/) 


0.178 



a Ephemeris taken from Kilkenny et al. (2000). 
b T c f{2 was kept fixed as it is poorly constrained by the data, see the 
text for details. 



both for the primary g[ and the secondary g2, bolometric albedo 
of the primary A 1 and the limb darkening coefficients of the pri- 
mary in the two filters x\ (g\ r). For the gravity darkening co- 
efficients we adopted values of 1.0 for the primary (radiative en- 
velope) and 0.32 for the secondary (convective envelope). We 
assumed a circular orbit (e=0) and synchronized rotation with 
the orbit. 

The effective temperature of the primary r e ffi was set to 
the value derived from our spectra (see Sect. |4.1[l. The effec- 
tive temperature previously estimated by Kilkenny et al. ( 1998, 
T e ff = 33 000 + 1 000) was used as well, but, as it did not in- 
fluence the derived parameters except for the luminosity of the 



M. Vuckovic et al.: The binary properties of the pulsating sdB eclipsing binary PG 1336-018 




HJD 



0.58 0.60 
2453509.0 



0.70 



Fig. 7. Same as Fig. [6] but for the f band. The trends seen in the middle and bottom panel result from imperfect removal of the 
reflection effect due to the changing temperature across the surface of the secondary (see text for details). 



stars, we fixed the temperature to the value derived by our new 
data. The T^n of the secondary has a very low contribution to the 
total flux (see Sect. |2]i and, therefore, is not tightly constrained. 
An appropriate treatment of the effective temperature of the sec- 
ondary in the case where the hot sdB primary is heating the cool 
secondary is not trivial, as the temperature on the illuminated 
hemisphere can be as much as five times higher than on the non- 
illuminated one dZolall2000l) . Whilst we did not intend to fix the 
effective temperature of the secondary star at first, we have found 
that leaving it as an adjustable parameter does not give consis- 
tent results. With T^i as a free parameter, it converges to around 
4000 K for the g' lightcurve, but to only 2700 K for the r'-band 
lightcurve. As a reasonable compromise for r e ff2, we choose to 
fix it to 3000 K. Considering the fact that the contribution of the 
secondary to the total flux is negligible, this is not an obstacle. 



As there are no published limb darkening coefficients for 
sdB stars we calculated the limb darkening coefficients x\ (g\ r* 
and m') for a 'typic al' sdB star from a fully line-blanketed LTE 
model atmosphere dBehara & Jeffery|2006l) with T eS = 30 000 K, 
logg = 5.5, V tul -b=5km/s and solar abundances (a linear cosine 
law was used). The mean limb darkening coefficients in each fil- 
ter were computed by convolving the ULTRACAM efficiencies in 
each filter with the monochromatic limb darkening coefficients 
and the stellar fluxes. We also computed the orbital solution 
using an ex trapolation of prev i ously repo rted coeffic i ents from 
the tables of lWade & Rucinskil d 1985b and lAl-Naimiyldl978l). as 
well a s the values fixed at 0.25 (V) and 0.20 (R) dKilkennv et al] 
1998). This did not change the solution, so we adopted the coef- 
ficients we computed from a modern atmosphere model. TableQ] 
summarises the values of the fixed parameters. The surface grav- 



ity is not a free parameter obtained by PHOEBE, since it is defined 
by the mass and radius. 



Using the ephemeris given in iKilkennv et ail d2000h we find 
a phase shift of 0.00374+0.00006 d. This phase shift could in 
principle be due to timing errors in our data rather than to an 
intrinsic change in the system. However, we carefully checked 
timings in our data sets and, moreover, we have data from two 
different instruments which both show the same phase shift. A 
timing error is therefore very unlikely to be the cause of the mea- 
sured shift. A change inherent to the system is thus the most 
probable reason. With only two minima timings we cannot draw 
any further conclusion here, only emphasise the need for fur- 
ther epoch observations. A similar period change on the order 
of 0.003 d over a period of 6 years i n the HW Vir system was 
documented by Kil kenny et al.l d2000h . 

The strong pulsations in the lightcurves are obstructing the 
fine tuning of the orbit, as the pulsations are seen as scatter by 
PHOEBE. Therefore, we take the first iteration solution and sub- 
tract it from the lightcurves. Now, after the dominant parts of 
the periodicity, i.e. the eclipses, have been removed from the 
lightcurves we can analyse them in order to take out the pul- 
sations of the primary from the lightcurves. 

A Fourier amplitude spectrum was calculated for each or- 
bit subtracted lightcurve to deduce the periodicities present in 
the data. The short timespan of our photometric data confines 
us with a frequency resolution of 54 fjHz. Since we are unable 
to resolve many of the closely spaced fre quencies in the spec- 
trum published by Kilk enny et alj d2003l) . we cannot use their 
peaks. We can only remove the periodicities we observe in our 
data in order to improve our orbit solution, after verifying that 
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Table 2. The list of frequencies, periods, amplitudes and phases we detected and prewhitened our data with. The phase is given as 
the time of maximum amplitude since to. 



Frequency 


Period 




Amplitude 






Phase (r raax ) 




[pHz] 


M 




[mma] 






[s] 








g' 


r 


m' 


g' 


r 


iC 


5430.1 


184.16 


11.2(1) 


10.5(1) 


17.1(2) 


142.3(3) 


142.2(3) 


141.4(4) 


5579.9 


179.21 


3.8(1) 


3.7(1) 


3.5(2) 


105.9(8) 


105.8(8) 


115(2) 


5757.3 


173.69 


1.7(1) 


1.7(1) 


2.8(2) 


148(2) 


147(2) 


155(2) 


7076.7 


141.31 


2.0(1) 


1.9(1) 


3.0(2) 


105(1) 


106(1) 


107(2) 



the frequencies we detect are indeed in the range of known 
PG 1336-018 frequencies. 

After identifying the highest amplitude peak in the spectrum 
and cross-checking if this frequency is present in the previous 
data sets within our frequency resolution, we remove this peak 
from the data by subtracting a sine wave (with the frequency, 
amplitude and phase determined by a non-linear least-squares fit 
-NLLS) from the original lightcurves. We calculate the Fourier 
amplitude spectrum of the prewhitened residuals and repeat the 
procedure until no new peaks could be securely identified. In this 
way we are able to remove four frequencies, as listed in Table[2] 
The frequency spectrum of PG 1336-018 is complicated as there 
are many frequencies in a narrow frequency range, which are 
unresolved in our data set. Therefore the NLLS would not con- 
verge on a simultaneous fit to more than four frequencies, even 
though there is still significant power left in the Fourier spec- 
trum. That is also the reason why the amplitude s appear higher 
in our data set compared to the ones seen in Kil kenny et alj 
(2003) as several frequencies are blended into one. The high- 
est amplitude frequency in our data set at 5430. 1 /zHz is most 
probably the result of seven unresolved closely s paced frequen- 
cies h , fj, fi5, j io, fs, fi and fn from Table 4 of iKilkennv et alj 
(12001 . 

These prewhitened lightcurves were then phase binned and, 
together with the RV curve, fed into PHOEBE to search for the 
improved orbit solution. Even though residual pulsations are still 
clearly visible in the lightcurves, their amplitudes are now signif- 
icantly smaller, which allows us to obtain a more reliable (sec- 
ond iteration) orbit solution. A third iteration step turns out to 
be unnecessary, as it does not improve the final outcome of the 
orbital parameters. 

As a quantitative measure of the goodness-of-fit we use the 
1 cr deviation for each data set (g\ r and RV) from the si- 
multaneously calculated synthetic curves. The bigger 1 cr de- 
viation in g' is due to the higher amplitudes of the oscillations 
in this colour. While it is impossible to see the depth of the lo- 
cal minima found by the DC method, and therefore search for 
the global minimum of the parameter hyperspace, we tested the 
stability of the conve rgent solutions found by parameter kicking 
dPrsa & Zwitterl2005l) . Once convergence was reached, we man- 
ually kicked the parameters and the minimization was restarted 
from the displaced points. In this way we found three groups of 
solutions of equal goodness-of-fit. Table [3] gives the three best 
fit orbital solutions. It is not possible to decide which solution is 
the correct one based on the numerical considerations as the syn- 
thetic curves are fitting the data equally well for all three models. 
The errors given in the table are the formal errors of the fit which 
are likely smaller than the true errors due to the above mentioned 
correlation between the parameters. The synthetic lightcurve fits 
to the observed data points are presented in Fig. [4] Fig. [6] and 
Fig- LZI (solid line) together with their residuals. The synthetic g' 
and f lightcurves and the RV curve are plotted for only one solu- 



Table 3. System parameters of the three best model fits to RV 
data and lightcurves of PG 1336-018. The formal lcr error on 
the last digit of each parameter is given in parentheses. 



Free parameter 


Model I 


Model II 


Model III 


a[R ] 


0.723(5) 


0.764(5) 


0.795(5) 


1 


0.282(2) 


0.262(2) 


0.250(2) 


hi 


80.67(8) 


80.67(8) 


80.67(8) 




5.50(3) 


5.48(3) 


5.47(3) 


n 2 


2.77(1) 


2.68(1) 


2.62(1) 


A 2 


0.92(3) 


0.92(3) 


0.93(3) 


*2 (g) 


0.38(8) 


0.39(8) 


0.38(8) 


x 2 (r') 


0.88(8) 


0.89(8) 


0.89(8) 


Derived parameters: 


Mi [Mo] 


0.389(5) 


0.466(6) 


0.530(7) 


M 2 [Mo] 


0.110(1) 


0.122(1) 


0.133(2) 


R, [Ro] 


0.14(1) 


0.15(1) 


0.15(1) 


R 2 [Ro] 


0.15(1) 


0.16(1) 


0.16(1) 


loggi [cm/s 2 ] 


5.74(5) 


5.77(6) 


5.79(7) 


logg 2 [cm/s 2 ] 


5.14(5) 


5.14(5) 


5.14(5) 


Roche radii: [in units of orbital separation] 


r, (pole) 


0.191 


0.191 


0.191 


r, (point) 


0.193 


0.193 


0.193 


r x (side) 


0.192 


0.192 


0.192 


ri (back) 


0.193 


0.193 


0.193 


r 2 (pole) 


0.198 


0.197 


0.197 


r 2 (point) 


0.213 


0.215 


0.216 


r 2 (side) 


0.201 


0.201 


0.201 


r 2 (back) 


0.210 


0.211 


0.211 


Errors on residuals: 


cr(gl [mag] 


0.03055 


0.03054 


0.03057 


cr(r') [mag] 


0.01325 


0.01321 


0.01321 


tr(RV) [km/s] 


8.39 


8.39 


8.39 



0.0014 




0.0011 ' ' ' ' ' ' ' 1 

0.22 0.24 0.26 0.28 0.3 0.32 

1 

Fig. 8. Mass ratio q versus sigma, for the range of the possible q 
values. Sigma is the sum of the squares of the sigmas in the two 
considered filters (cr(g') and cr(r')). 
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Fig. 9. Mass-radius diagram for PG 1336-018 showing the re- 
gions permitted by the orbit solution (continuous line) and by the 
different surface gravities (dotted lines). The q values are also 
noted on the orbit solution. The small changes from the 3<x error 
on Ki do not shift the curve representing the orbital solution. 

tion (Model II) since the deviations between the three solutions 
cannot be resolved at the scale of the figure. 

4.3. Discussion 

The uniqueness of a given solution is jeopardized by the param- 
eter correlations. In particular, there is a strong correlation be- 
tween the mass ratio q and the potential of the secondary star 
£^2- Therefore, there is a q degeneracy in all the orbital solutions. 
For a given range of potentials defined by the Lagrangian point, a 
family of solutions with corresponding mass ratios is found. The 
solutions found in Table [3] represent the local minima shown in 
Fig.E] 

The relative radii and the orbital inclination are tightly con- 
strained by the depth and the width of the eclipses, and the re- 
sults in all three models are nearly identical. There is only a 
slight distortion of the secondary: r2 (pole)/r2 (point) is 0.93, 
0.92, 0.91 respectively for each model . While the previous 
searches for the be st orbital solutions dKilkennv et al.l [1998; 
iDrechsel et al.l l200ll and references therein) tend to resort to 
non-physical albedos (greater that 1 in some cases) and limb 
darkening coefficients of the secondary, we find that the biggest 
problem is in the temperature of the secondary which is heated 
by the hot subdwarf. The weakest point of all modelling proce- 
dures lies in an inadequate treatment of the temperature of the 
secondary star. The temperature distribution over the surface of 
the secondary has to be incorporated in the atmosphere models 
used by PHOEBE in order to get more realistic solutions. This is 
far beyond the scope of our current paper. 

The surface gravity derived from the orbital solutions, al- 
though in agreement with the value previously estimated by 
Kilk enny et all dl998l logg = 5.7 ± 0.1 dex) is higher than the 
spectroscopic gravity estimate. Therefore, we have explored the 
full range of mass-radius ranges for the primary allowed by the 
orbital solution and the spectroscopic gravity (Fig. [9]). The pa- 
rameters used to generate this orbital solution mass-radius re- 
lationship are only the P, i, K\ and the radius of the primary 
in terms a, none of which are affected by the q degeneracy. 
Thus, if we had a sufficiently accurate spectroscopic determi- 



25 



20 - 




-20 - 



-25 1 ' ' 1 ' 1 ' ' 1 

-0.15 -0.1 -0.05 0.05 0.1 0.15 

Phase 

Fig. 10. The orbit subtracted RV residuals (dots) with their cor- 
responding errors clearly showing the RM effect. The solid line 
is the simulation of the RM effect with the parameters given in 
the text. 

nation of log g, we could use the relationships in Fig.|9]to deter- 
mine one unique M\. Unfortunately, our spectroscopic \ogg of 
5.6 is clearly much lower than what can realistically be accepted 
since it gives a mass for the primary that is far too low (M\ < 0.2 
[Mo]). 

While we cannot discriminate between the three model fits 
on the basis of their <x values, the evolutionary scenarios for sdB 
stars disqualify the Model III solut ion as the prim ary mass would 
be too high for a core He-flash (lHan et al.| [2002). Models I and 
II however, are both possible as they could have formed through 
common envelope phase (Hu et al., submitted to A&A). 

5. Detection of the Rossiter-McLaughlin effect 

In Fig. [4] an apparent up-and-down (redshift-blueshift) shift oc- 
curs at phase zero in the RV curve. This effect at the eclipse is 
known as the Ross iter-McLaughlin (RM) effect dRossiterlll924t 
iMcLaughlinl 119241) . It is due to the selective blocking of the 
light of the rotating star during an eclipse. When the secondary 
star covers the blueshifted (redshifted) half of the stellar disk, 
the integrated light of the primary appears slightly redshifted 
(blueshifted). Because of this selective blocking of the stellar 
surface during the eclipse, a skewed line profile is created. This 
change in line profile shape results in a shift in RV, which in turn 
results in the redshift-blueshift distortion seen during the eclipse 
(see Fig. 3J. The RM effect has bee n seen in other eclipsin g 
hot subdwarf binaries (e.g. AA Dor: iRauch & Wernerl (120031) ) 
and can be used to investigate the rotational properties of the 
compon ent stars. It was recently used in ext r asolar planetary 
transits dOueloz et all 120001 lOhta et al.l |200H iGimenezI 120061: 
Gaudi & Winn 120061) to discriminate between different migra- 
tion theories. The amplitude of the effect mainly depends on the 
projected rotation velocity of the star, the ratio of stellar radii, 
the orbital inclination, and the limb darkening. 

To analyze this effect we have subtracted the orbital solution 
(solid curve in Fig. |4]i from the RV measurements. The orbit- 
subtracted RV residuals, phase binned in 50 bins, are plotted in 
Fig. [10] The RM effect is clearly seen in these r esiduals. We used 
the analytical description of this effect given in Gimenez (2006) 
to simulate the RM effect for this system. We have assumed that 
the rotational axis of the primary star is co-aligned with the per- 
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pendicular to the orbital plane. The result of this simulation is 
plotted as a solid line in Fig. [10] The equatorial rotational veloc- 
ity of the star was set to 75.2 kms -1 and the ratio of the stellar 
radii rijrx , the inclination of the orbit i and the radius of the pri- 
mary relative to the size of the orbit r\ were taken from our or- 
bital solution (see Table |3). The synthetic curve fits the observed 
RM amplitude rather well. The uncertainties on the residual RV 
curve are too large to fine-tune the orbital parameters. We can 
only establish that the observed RM effect is compatible with 
the orbital solutions given in Tableland represents an indepen- 
dent confirmation of the light curve solution. 

The apparent asymmetry seen in Fig.[l0]is, however, not well 
explained. Such an asymmetry is expected to occur if the pro- 
jected orbital and rotational axes are not aligned. This is highly 
unlikely for the narrow orbit of PG 1336-018. Nevertheless, we 
simulated the RM effect allowing different angles of the rotation 
axes and the orbital axes. We indeed could not achieve satisfac- 
tory results, because, when the zero offset was fitted well, the 
amplitudes were highly asymmetrical and vice versa. The asym- 
metry is more likely caused by the pulsations seen during the 
primary eclipse, which also give rise in line profile shape vari- 
ations. The equations describing the RM effect assume that the 
components are spherical, i.e. they do not take into account any 
deviation from spherical symmetry such as the one produced by 
the pulsations. We will investigate this further in our follow-up 
paper dedicated to the analysis of the primary's pulsations. 

6. Conclusions and Future work 

In this work, we presented a thorough observational analy- 
sis of the orbital behavior of the pulsating eclipsing binary 
PG 1336-018. Our goal was to avoid using a canonical mass 
of 0.5 M for the subdwarf in any interpretation of the luminos- 
ity variations of the star, as has been done so far in the litera- 
ture. Instead, we attempted an unbiased derivation of the system 
and stellar parameters, in particular for the masses of the com- 
ponents. Our analysis resulted in three equally probable sets of 
orbital and physical parameters of the system. Our model III so- 
lution is incompatible with the binary having gone through a 
core He-flash and a common-envelope phase described by the 
a-formalism since that can only lead to PG 1336-018 like bina- 
ries with primary masses up to 0.48 M (Hu et al., submitted to 
A&A). This leaves us with two solutions, one with a primary 
mass of 0.466+0.006 M and another with 0.389+0.005 M , 
with secondary masses of 0. 1 22+0.00 1 M and 0. 1 1 0+0.00 1 M 
respectively. We thus conclude that our solutions with M\ = 
0.466+0.006 M and M x = 0.389+0.005 M are the only plau- 
sible ones, except when the common- envelope phase would 
be better described by the y-formalism (Nele mans et al.|[2000t 
iNelemans & Tout! [20051) . In this case all three solutions are ac- 
ceptable, as this formalism allows non-degenerate helium igni- 
tion with a broader primary mass range (0.3-1.1 M ). 

Furthermore, we have detected the RM effect in the radial 
velocity curve of PG 1336-018. The simulated amplitude of the 
RM effect is in the accordance with the RM amplitude seen in 
the RV residuals, which is an independent confirmation of the 
results obtained from our orbital solution. 

While deriving the orbital solution for PG 1336-018, we hit 
upon the limitation of current binary analysis codes, which also 
prevented us to pinpoint the effective temperature of the sec- 
ondary. None of the analysis methods available in the literature 
treat the atmosphere of such a close binary, in which one com- 
ponent is so hot that it induces a temperature gradient across 
the surface of the other, in an appropriate way. Indeed, all codes 



make use of stellar atmosphere models which assume one fixed 
effective temperature at the surface of each of the component 
stars. As such, any derived quantities, such as limb darkening 
coefficients and albedos, cannot be but a very crude approxima- 
tion of reality whenever one component is seriously heated by 
the other one. In the case of close binaries like PG 1336-018, 
i.e. with a hot primary and a cold secondary, the temperature 
of the latter changes so drastically from the illuminated side to 
the backside, that specific atmosphere models representing such 
a situation should be computed and used while deriving the or- 
bital parameters. This is an entire project by itself and surely 
beyond the scope of our current work. We hope that our results 
will give rise to future developments of atmosphere models with 
temperatures varying across the surface of the cool component 
in close binaries. The case of PG 1336-018, and our data of the 
star, are ideally suited to test such new future models. 

In a follow-up paper of this work, we plan to analyse 
the oscillatory signal in our multicolour photometry and high- 
resolution spectroscopy, after the orbit subtraction presented 
here. This will be done by computing a cross-correlation func- 
tion of each spectrum and investigating the signature of the 
modes in it. Cross-correlation functions have already been 
use d to study the character of oscillations modes before, see 
e.g iMathias & Aertsl d 19961) for the 5Scuti star 20CVn and 
Hek ker et al.l ([2006) for solar-like oscillations in red giants. This 
is done by computing line diagnostics, such as moments, and 
the amplitude and phase across the profile, and comparing these 
to predictions based on the theory of non-radial oscillations. In 
principle, this allows us to identify the spherical wavenumbers 
{(, m) of the strongest modes. The use of these established mode 
identi fication techniques (see e.g. iBriquet & Aertsl 12003b IZimal 
2006, for the latest versions) on high-resolution cross-correlation 
profiles of pulsating sdB stars has so far not yet been done. The 
nature of our data and of our target star requires a simulation 
study to test the effects of smearing out the oscillations over 
the cycle and of the limited time base. Also, we must treat the 
data during and outside the eclipses separately in order to as- 
sess the effectiveness of the techniques in the specific case of 
PG 1336-018. Such a study is currently being performed. The 
ultimate goal of it is to identify the highest-amplitude modes 
and discriminate among the plausible seismic models of the star. 
This will then eventually lead us to derive a seismic mass es- 
timate to be confronted with the observed primary masses pre- 
sented here and with the evolutionary masses computed by Hu 
et al. (submitted to A&A). 
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